function [HH,stdVec,flagNotPD]=...
    modeFinderHessian(gSteps,parStru,posStru,model,dataStru,prpar,prss,filterStru,flags)
% Compute the Hessian through Finite Point Difference
% Called from modeFinderOutput.m 
%% model 
%  .handle 
%  .addsol 
%  .solveOptions  
%  .param           complete parameter vector 
%  .logPost         log Posterior at model.param
%% All remaining entries are to be used with modePosterior.m 
% see that function for details 
% A. Justiniano 
if isempty(gSteps)==true
    gSteps=1e-5;
end
modeEstimatedOnly=model.param(posStru.param.est); 
dim.est=length(modeEstimatedOnly);
disp('Begin Hessian');
%% 1. Check Initial Density Matches Posterior at Mode 
fInitial=...
    feval(filterStru.funcPost,modeEstimatedOnly,parStru,...
    posStru,model,dataStru,prpar,prss,filterStru,flags);
fprintf('Starting Log Posterior for Hessian=%10.5f \n',fInitial);
if abs(fInitial-model.logPost) > 1e-4
    warning('Starting Posterior for Hessian and Posterior at Mode do not coincide')
end
%% 2. Beging Hessian Routine 
HH=hessian(filterStru.funcPost,modeEstimatedOnly,gSteps,...
    parStru,posStru,model,dataStru,prpar,prss,filterStru,flags);
disp('Done with Hessian')
try
    save Hessian
catch
    disp('Could not save Hessian')
end
%% 3. Attempt Inversion by SVD 
HH=reshape(HH,[dim.est dim.est]);
HH=generalized_cholesky(HH);
% AJ Feb 2014 
HH=svdInverse(HH,1e-10); 
[~,flagNotPD]=chol(HH);
if flagNotPD~=0
    warning('Hessian is not positive definite')
    flagNotPD=1;
    stdVec=nan(dim.est,1); 
    return 
end
stdVec=sqrt(diag(HH));